# 2023.9.15
library(Seurat)
library(tidyseurat)
library(tidySingleCellExperiment)
library(ggpubr)
source("00_util_scripts/mod_seurat.R")

reduce <- purrr::reduce
batch.c4 <- c('orig.ident','dataset')

# load seurat objects ----------
sobj_sg <- read_rds('CRC-I/data/seekgene/crc2305.crEM.rds')

sobj_10x <- read_rds('CRC-I/data/zy2020_tumor10x.rds')

sobj_guo <- read_rds("CRC-I/data/guo2021sided/guo2021_immune.rds")

sobj_liv <- read_rds('CRC-I/data/li_2021_liver/li2021crc.rds')

# merge 4 datasets ----------
sobj_guo$dataset <- 'guo2021'
sobj_liv$dataset <- 'li2021'
sobj_sg$dataset <- 'liu2023'
sobj_10x$dataset <- 'zhang2020'

sobj_liv %<>%
  mutate(genotype = FCGR2B.I232T)

# remove non-immune cells
sobj_liv <- sobj_liv |>
  filter(!str_detect(hpca_main,'Epi|stem|Endo'))

gene_merged <- c(sobj_sg,sobj_liv,sobj_guo,sobj_10x) |>
  map(rownames) |>
  reduce(intersect)

# this can take ~1 min
sobj_c4 <- c(sobj_sg, sobj_10x, sobj_liv, sobj_guo) |>
  reduce(merge) |>
  subset(features = gene_merged) |>
  JoinLayers()

sobj_c4 <- sobj_c4 |>
  select(.cell, orig.ident:nFeature_RNA,
         genotype, dataset)

# 85k cells, take ~5.3 min for 25 pcs
# ~12.6 min for 100 pcs
sobj_c4 <- sobj_c4 |>
  quick_process_seurat(batch = c('orig.ident','dataset'),
                       pcs = 25,
                       skip_norm = TRUE)

sobj_c4 |> DimPlot(split.by = 'dataset',ncol = 2)

## manually annotate ----------
sobj_c4 |>
  DotPlot(features = seurat_markers, cluster.idents = T) +
  RotatedAxis()

### marker for ambiguous clusters -------
ambi_markers <- c(11,16:21) |>
  map(\(x)FindMarkers(sobj_c4, ident.1 = x, only.pos = TRUE, logfc.threshold = 2)) |>
  map(rownames_to_column, 'gene') |>
  set_names(c(11,16:21)) |>
  list_rbind(names_to = 'cluster')

ambi_markers |>
  filter(cluster == 16) |>
  slice_max(avg_log2FC, n = 10)

no_other_b_11 <- sobj_c4 |>
  filter(seurat_clusters != 3 & seurat_clusters != 4) |>
  FindMarkers(ident.1 = 11, only.pos = TRUE, logfc.threshold = 1.5)

# comparing dotplot and dimplot to manual annotate merge4, do not delete this!
sobj_c4 <- sobj_c4 |> mutate(manual_main = case_match(
  as.numeric(seurat_clusters),
  c(2) ~ 'B',
  c(4,14) ~ 'Plasma',
  c(1,5,10,12) ~ 'CD4T',
  c(2,6,18) ~ 'CD8T',
  c(7,9,14:17) ~ 'Myeloid',
  9 ~ 'Mast',
  12 ~ 'NK',
  11 ~ 'Dividing lymphocyte',
  c(19,21) ~ 'Epithelia'))

sobj_c4 <- sobj_c4 |>
  filter(!is.na(manual_main))

sobj_c4 |> DimPlot(group.by = 'manual_main', label = TRUE,
                   label.box = TRUE) +
  ggtitle('Major immune cells')

sobj_c4 |> DimPlot(group.by = 'manual_main', split.by = 'dataset',ncol = 2) +
  ggtitle('Major immune cells')

## auto annotate ------
sobj_c4 |> DimPlot(group.by = 'hpca_main', label = TRUE,
                   label.box = TRUE)

monaco <- celldex::MonacoImmuneData()

sobj_c4 |>
  filter(seurat_clusters %in% c(3,4,11)) |>
  mark_cell_type_singler(monaco, fine_label = TRUE, new_label = 'monaco_fine')

hpca <- celldex::HumanPrimaryCellAtlasData()

sobj_c4 <- sobj_c4 |>
  mark_cell_type_singler(hpca, new_label = 'hpca_main_c4')

sobj_c4 |>
  DimPlot(group.by = 'hpca_main_c4', cols = 'Paired')

crc10x.ref <- read_rds('CRC-I/data/crc10x_singler_ref.rds')

sobj_c4 <- sobj_c4 |>
  mark_cell_type_singler(crc10x.ref, new_label = 'zhang.main')

sobj_c4 |>
  DimPlot(group.by = 'zhang.main', cols = 'Paired')

sobj_c4 <- sobj_c4 |>
  mark_cell_type_singler(crc10x.ref, new_label = 'zhang.fine',
                         sc_ref = TRUE, fine_label = T)

sobj_c4 |>
  DimPlot(group.by = 'zhang.fine', cols = DiscretePalette(36))

## final type -------
sobj_c4 %<>%
  mutate(manual.main = case_when(
    seurat_clusters == 9 ~ 'Mast_cells',
    hpca_main_c4 %in% c('T_cells','NK_cell') ~ hpca_main_c4,
    seurat_clusters == 2 ~ 'B_cells',
    seurat_clusters %in% c(4,15) ~ 'Plasma_cells',
    .default = 'Myeloid_cells'
  ))

sobj_c4 |> DimPlot(group.by = 'manual.main',
                   cols = DiscretePalette(36))+
  ggtitle('Major immune cells')

sobj_c4 |> DimPlot(group.by = 'hpca_main_c4', split.by = 'dataset',ncol = 2,
                   cols = DiscretePalette(10)) +
  ggtitle('Major immune cells')

### eliminate ambiguous & non-immune cells -------
marker11 <- FindMarkers(sobj_c4, ident.1 = 11, only.pos = TRUE, logfc.threshold = 1)

marker11 |>
  as_tibble(rownames = 'gene') |>
  arrange(desc(avg_log2FC))

# outstanding T cell cluster are cycling T! ----
DotPlot(sobj_c4, features = cc.genes.updated.2019)

sobj_c4 |> write_rds('CRC-I/data/crc_merge4_immune.rds')

sobj_c4 <- read_rds('CRC-I/data/crc_merge4_immune.rds')

sobj_c4@meta.data |>
  as_tibble(rownames = '.cell') |>
  write_csv('CRC-I/results/crc_merge4_meta.csv')

mast_tt_deg <- sobj_c4 |>
  filter(manual_main == 'Mast') |>
  FindMarkers(ident.1 = 'TT', group.by = 'genotype')

nk_tt_deg <- sobj_c4 |>
  filter(manual_main == 'NK') |>
  FindMarkers(ident.1 = 'TT', group.by = 'genotype')

main_tt_deg <-sobj_c4$manual_main |>
  unique() |>
  map(\(x)sobj_c4 |> filter(manual_main == x) |>
        FindMarkers(ident.1 = 'TT', group.by = 'genotype')) |>
  map(rownames_to_column, 'gene') |>
  set_names(unique(sobj_c4$manual_main)) |>
  list_rbind(names_to = 'cluster')

main_tt_deg |> write_csv('CRC-I/results/merge4_main_tt_deg.csv')

# subclusters ------
## myeloid --------
sobj.myl <- sobj_c4 |>
  filter(str_detect(manual.main, 'Myel|Mast'))

# cost 67s
sobj.myl %<>%
  quick_process_seurat(batch = batch.c4,
                       skip_norm = T)

sobj.myl %<>%
  mark_cell_type_singler(
    ref = filter(crc10x.ref, str_detect(label.main,'Myel')),
    sc_ref = TRUE,
    fine_label = TRUE,
    new_label = "zhang.fine"
  )

sobj.myl |> DimPlot(group.by = 'zhang.fine', cols = 'Paired') +
  ggtitle('CRC TME myeloid cells')

sobj.myl %<>%
  mutate(namely.fine = case_when(
    seurat_clusters == 10 ~ 'Monocyte',
    str_detect(zhang.fine, 'Mast') ~ 'Mast cell',
    str_detect(zhang.fine, 'DC') ~ 'DC',
    .default = 'Macrophage'
  ))

sobj.myl |> DimPlot(group.by = 'namely.fine') +
  ggtitle('CRC TME myeloid cells')

sobj.myl |>
  as_tibble() |>
  select(.cell, zhang.fine, namely.fine) |>
  write_csv('CRC-I/results/crc_merge4_myel.meta.csv')

## B & PC ---------
sobj.bc <- sobj_c4 |>
  filter(str_detect(manual.main, 'B|Plasma'))

# cost 59s
sobj.bc %<>%
  quick_process_seurat(batch = batch.c4,
                       skip_norm = T)

sobj.bc %<>%
  mark_cell_type_singler(
    ref = crc10x.ref,
    sc_ref = TRUE,
    fine_label = TRUE,
    new_label = "zhang.fine"
  )

sobj.bc %<>%
  filter(seurat_clusters != 13)

sobj.bc |> DimPlot(group.by = 'zhang.fine', cols = 'Paired') +
  ggtitle('CRC TME B cells')

sobj.bc %<>%
  mutate(namely.fine = case_when(
    str_detect(zhang.fine, 'Follicu') ~ 'B cell',
    str_detect(zhang.fine, 'GCB') ~ 'GC B cell',
    .default = 'Plasma Cell'
  ))

sobj.bc |> DimPlot(group.by = 'namely.fine') +
  ggtitle('CRC TME B cells')

sobj.bc |>
  as_tibble() |>
  select(.cell, zhang.fine, namely.fine) |>
  write_csv('CRC-I/results/crc_merge4_bc.meta.csv')

## T+NK -------
sobj.nkt <- sobj_c4 |>
  filter(str_detect(manual.main, 'T|NK'))

# cost 59s
sobj.nkt %<>%
  quick_process_seurat(batch = batch.c4,
                       skip_norm = T)

sobj.nkt %<>%
  mark_cell_type_singler(
    ref = crc10x.ref,
    sc_ref = TRUE,
    fine_label = TRUE,
    new_label = "zhang.fine"
  )

sobj.nkt %<>%
  filter(!str_detect(zhang.fine, 'BCell|Plasma'))

sobj.nkt |> DimPlot(group.by = 'zhang.fine', cols = 'Paired') +
  ggtitle('CRC TME T/NK cells')

sobj.nkt %<>%
  mutate(namely.fine = case_when(
    str_detect(zhang.fine, 'CD8') ~ 'CD8 T cell',
    str_detect(zhang.fine, 'NK') ~ 'NK cell',
    str_detect(zhang.fine, 'ILC3') ~ 'ILC3',
    .default = 'CD4 T Cell'
  ))

sobj.nkt |> DimPlot(group.by = 'namely.fine') +
  ggtitle('CRC TME T cells')

sobj.nkt |>
  as_tibble() |>
  select(.cell, zhang.fine, namely.fine) |>
  write_csv('CRC-I/results/crc_merge4_nkt.meta.csv')

## examine FCGR2B expression ---------
fcgr2b_expr <- sobj_c4 |>
  get_abundance_sc_long('FCGR2B')

fcgr2b_expr |>
  filter(.abundance_RNA > 0) |>
  dplyr::select(.cell) |>
  write_csv('CRC-I/results/merge4.fcgr2b.pos.cell.csv')

fcgr2b_expr <- sobj_c4 |> as_tibble() |>
  dplyr::select(.cell, manual.main) |>
  left_join(fcgr2b_expr)

fcgr2b_expr |>
  ggplot(aes(manual.main, .abundance_RNA, color = manual.main)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white', linewidth = 2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_cl_boot', width = .3) +
  ggpubr::theme_pubr() +
  labs(x = 'Cell type', y = 'Normalized expression', color = 'Cell type',
       title = 'FCGR2B expression') +
  NoLegend()

g1 <- last_plot()

fcgr2b_expr |>
  group_by(manual.main) |>
  mutate(wt = 1/n()) |>
  filter(.abundance_RNA > 0) |>
  dplyr::count(manual.main, wt = wt) |>
  ggplot(aes(manual.main, n, fill = manual.main)) +
  geom_col() +
  ggpubr::theme_pubr() +
  labs(x = 'Cell type', y = 'Fraction of FCGR2B+ cells', color = 'Cell type',
       title = 'Fraction of FCGR2B+ cells') +
  NoLegend()

g2 <- last_plot()

g1 / g2 + patchwork::plot_annotation(title = 'CRC TME immune cell expression of FCGR2B')

### only cd8 -----
fcgr2b_expr <- sobj.cd8 |>
  get_abundance_sc_long('FCGR2B')

fcgr2b_expr <- sobj.cd8 |> as_tibble() |>
  select(.cell, zhang2020_main) |>
  left_join(fcgr2b_expr)

fcgr2b_expr |>
  ggplot(aes(zhang2020_main, .abundance_RNA, color = zhang2020_main)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white', linewidth = 2) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_cl_boot', width = .3) +
  ggpubr::theme_pubr() +
  labs(x = 'Cell type', y = 'Normalized expression', color = 'Cell type') +
  NoLegend() +
  coord_flip() +
  scale_color_brewer(type = 'qual')

g1 <- last_plot()

fcgr2b_expr |>
  group_by(zhang2020_main) |>
  mutate(wt = 1/n()) |>
  filter(.abundance_RNA > 0) |>
  count(zhang2020_main, wt = wt) |>
  ungroup() |>
  add_case(zhang2020_main = 'LAYN-CD8-Tex', n = 0) |>
  ggplot(aes(zhang2020_main, n, fill = zhang2020_main)) +
  geom_col() +
  ggpubr::theme_pubr() +
  labs(x = 'Cell type', y = 'Fraction of FCGR2B+ cells', color = 'Cell type') +
  NoLegend()+
  coord_flip() +
  scale_fill_brewer(type = 'qual')

g2 <- last_plot()

g1 + g2 + patchwork::plot_annotation(title = 'CRC TME CD8 T cell expression of FCGR2B (Guo, 2021, JCI Insight)')

sobj |> write_rds("CRC-I/data/guo2021sided/guo2021_ident.rds")

# VEGFA/TNF+ Mast cell -------------
sobj.c4.mast <- sobj_c4 |>
  filter(hpca_main_c4 == 'Mast cell')

sobj.c4.mast %<>%
  quick_process_seurat(batch = c('orig.ident','dataset'),
                       skip_norm = T)

sobj.c4.mast |>
  DotPlot(c('VEGFA','TNF'))

mast.dict <- sobj.c4.mast |>
  FindAllMarkers(features = c('VEGFA','TNF'), only.pos = T)

mast.dict

sobj.c4.mast %<>%
  mutate(mast.subtype = case_when(seurat_clusters %in% c(2,3) ~ 'VEGFA+',
                                  seurat_clusters %in% c(12,10) ~ 'TNF+',
                                  .default = 'Other'))

sobj.c4.mast |>
  DimPlot(group.by = 'mast.subtype')

sobj.c4.mast@meta.data |>
  as_tibble(rownames = '.cell') |>
  write_csv('CRC-I/results/crc_merge4_mast_meta.csv')
